rm(list=ls())

library(rstan) 
library(ggplot2)
library(Hmisc)
library(foreign)
library(car)
library(lme4) 
library(shinystan)
library(readr)

multi_final<-  read_csv("(...)/pseudo_final.csv")

#View(multi_final)
attach(multi_final)


# Create a non-missing dataset for education

subset_health <- c("cnt", "av_spend_health","std_av_party", "std_av_income", "std_gini_net", "std_unemp", "std_gdp_growth","std_laghealth_av","groupid", "year")
multi_final_health <- multi_final[subset_health]    
#View(multi_final_edu)



multi_final_health2 <-   na.omit(multi_final_health)
#View(multi_final_health2)
attach(multi_final_health2)


#Prepare the variables for the analysis

recyear <- as.numeric(year == 2006)
y <- multi_final_health2$av_spend_health
summary(y)
groupid1 <- as.numeric(groupid == 1)
groupid2 <- as.numeric(groupid == 2)
groupid3 <- as.numeric(groupid == 3)
groupid4 <- as.numeric(groupid == 4)
groupid5 <- as.numeric(groupid == 5)
groupid6 <- as.numeric(groupid == 6)
groupid7 <- as.numeric(groupid == 7)
groupid8 <- as.numeric(groupid == 8)
groupid9 <- as.numeric(groupid == 9)
groupid10 <- as.numeric(groupid == 10)
groupid11 <- as.numeric(groupid == 11)
groupid12 <- as.numeric(groupid == 12)
groupid13 <- as.numeric(groupid == 13)
groupid14 <- as.numeric(groupid == 14)
groupid15 <- as.numeric(groupid == 15)
groupid16 <- as.numeric(groupid == 16)
groupid17 <- as.numeric(groupid == 17)
groupid18 <- as.numeric(groupid == 18)
groupid19 <- as.numeric(groupid == 19)
groupid20 <- as.numeric(groupid == 20)
groupid21 <- as.numeric(groupid == 21)
groupid22 <- as.numeric(groupid == 22)
groupid23 <- as.numeric(groupid == 23)
groupid24 <- as.numeric(groupid == 24)
groupid25 <- as.numeric(groupid == 25)
groupid26 <- as.numeric(groupid == 26)
groupid27 <- as.numeric(groupid == 27)

cnt2<- as.numeric(cnt)
party <- round(std_av_party,digits=2) 
income <- round(std_av_income,digits=2)
gini <- round(std_gini_net,digits=2)
unemployment <- round(std_unemp,digits=2)
laghealthav <- round(std_laghealth_av,digits=2)
growth <- round(std_gdp_growth,digits=2)
J<-21
attach(multi_final_health2)


# 1st part of the Bayesian model, define the data
modelhealthav = "
data {
int<lower=0> N; 
int<lower=1> J; 
int<lower=1,upper=J> cnt[N];
real <lower=-4.35,upper=4.44> party[N];
real <lower= -2.46,upper=2.81 >income[N];
real <lower=-1.59,upper=2.02> gini[N];
real <lower=-1.29,upper=2.52 > unemployment[N];
real<lower=-1.87,upper=3.04 > growth[N];
real<lower=-1.70,upper=2.41> laghealthav[N];
int<lower=0,upper=1> recyear[N];
int<lower=0,upper=1> groupid2[N];
int<lower=0,upper=1> groupid3[N];
int<lower=0,upper=1> groupid4[N];
int<lower=0,upper=1> groupid5[N];
int<lower=0,upper=1> groupid6[N];
int<lower=0,upper=1> groupid7[N];
int<lower=0,upper=1> groupid8[N];
int<lower=0,upper=1> groupid9[N];
int<lower=0,upper=1> groupid10[N];
int<lower=0,upper=1> groupid11[N];
int<lower=0,upper=1> groupid12[N];
int<lower=0,upper=1> groupid13[N];
int<lower=0,upper=1> groupid14[N];
int<lower=0,upper=1> groupid15[N];
int<lower=0,upper=1> groupid16[N];
int<lower=0,upper=1> groupid17[N];
int<lower=0,upper=1> groupid18[N];
int<lower=0,upper=1> groupid19[N];
int<lower=0,upper=1> groupid20[N];
int<lower=0,upper=1> groupid21[N];
int<lower=0,upper=1> groupid22[N];
int<lower=0,upper=1> groupid23[N];
int<lower=0,upper=1> groupid24[N];
int<lower=0,upper=1> groupid25[N];
int<lower=0,upper=1> groupid26[N];
int<lower=0,upper=1> groupid27[N];
vector[N] y;
} 

parameters {
vector[J] a;
vector[33] beta;
real mu_a;
real<lower=0,upper=10> sigma_a;
real<lower=0,upper=10> sigma_y;

}



transformed parameters {

vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- a[cnt[i]] +beta[1]*party[i]+ beta[2]*income[i] + beta[3]*gini[i]+ beta[4]*unemployment[i] + beta[5]*growth[i]+beta[6]*laghealthav[i]+beta[7]*recyear[i]+beta[8]*groupid2[i]+beta[9]*groupid3[i]+
beta[10]*groupid4[i]+beta[11]*groupid5[i]+beta[12]*groupid6[i]+beta[13]*groupid7[i]+beta[14]*groupid8[i]+beta[15]*groupid9[i]+
beta[16]*groupid10[i]+beta[17]*groupid11[i]+beta[18]*groupid12[i]+
beta[19]*groupid13[i]+beta[20]*groupid14[i]+beta[21]*groupid15[i]+beta[22]*groupid16[i]+beta[23]*groupid17[i]+beta[24]*groupid18[i]+beta[25]*groupid19[i]+beta[26]*groupid20[i]+
beta[27]*groupid21[i]+beta[28]*groupid22[i]+beta[29]*groupid23[i]+beta[30]*groupid24[i]+beta[31]*groupid25[i]+beta[32]*groupid26[i]+beta[33]*groupid27[i];
} 

model {
mu_a ~ normal(0, 10);
a ~ normal(mu_a, sigma_a);
beta ~ normal(0, 10);
y ~ normal(y_hat, sigma_y);
sigma_a ~ cauchy(0, 2.5);
sigma_y ~ cauchy(0, 2.5);
}

"
sm_healthav <- stan_model(model_code= modelhealthav)


dataList.1 <- list(N=length(y),y=y,J=J, cnt=cnt,party=party,income=income,gini=gini, unemployment=unemployment,growth=growth,laghealthav=laghealthav,groupid=groupid, recyear=recyear)

fit_healthav <- sampling(sm_healthav, data=dataList.1, iter=5000,chains=5)


shinystan_healthav <- as.shinystan(fit_healthav)
launch_shinystan(shinystan_healthav)



